!!         eval_jac.f90
!! This subroutine evaluates the Jacobean in the omega-psi
!! formulation of the 2-d navier-stokes equation. 
!! This is a serial code.This IS the part that contains
!! call to fourier transform subroutines. This SHOULD be changed
!! from one serial machine to another.
!! -----------------------------------------------------------------
subroutine eval_jac(flag,idiag) 
  use general
  use mod_2dflu
  use mod_part_interp   
  use fft
  use force
  use fourier_space_diagnostic
  use immersed_boundary
  implicit none
  double precision :: re_omega,im_omega,jplus,jminus,enx,eny,enkx,enky,ekx,eky,sum
  double precision, dimension(3) :: uu=0.,ww=0.,wXu=0.
  integer ::i1,i2
  integer:: flag,idiag
!! --------------------------------------------------------------
!! We recognise that the two derivative of psi appearing in the
!! Jacobean are simply the velocities. Hence we first evaluate
!! the velocities in Fourier space
  call calc_fspace_quant(idiag,flag)
!! ----------------------------------------------------
!! Transform velocities and vorticity to real space.
!!-----------------------------------------------------
  call fft_2d_complex2real(ukx)
  call fft_2d_complex2real(uky)
  call fft_2d_complex2real(omega)
  do i1=n1+1,n1+2
    ukx(i1,:) = 0.0d0
    uky(i1,:) = 0.0d0
    omega(i1,:) = 0.0d0
  enddo
  ukx = ukx*scale
  uky = uky*scale
  omega = omega*scale
!
!! ------------The velocities are in real space. ------------------
!
!!
!! Also get Gradu to Real space 
!!
  if((lcalc_Gradu).and.(idiag.eq.1).and.(flag.eq.1)) then
    call fft_2d_complex2real(dx_ux)
    call fft_2d_complex2real(dy_ux)
    call fft_2d_complex2real(dx_uy)
!!
    do i1=n1+1,n1+2
      dx_ux(i1,:) = 0.0d0
      dy_ux(i1,:) = 0.0d0
      dx_uy(i1,:) = 0.0d0
    enddo
    dx_ux = dx_ux*scale
    dy_ux = dy_ux*scale
    dx_uy = dx_uy*scale
!!
  endif
!
! calculate real space diagnostic here (including real space structure function)
! Also calculate the rhs of the particle equation and the immersed boundary eqn. 
!
  if (flag.eq.1) then 
    if (lparticle) call get_particle_rhs(ukx(1:n1,1:n2),uky(1:n1,1:n2),n1,n2)
    if (limmersed_boundary) call get_immersed_boundary_rhs(ukx(1:n1,1:n2),uky(1:n1,1:n2),n1,n2)
    if((idiag.eq.1)) then
!      call force_inj
      if(lparticle) call lambda_at_particle_position(dx_ux(1:n1,1:n2),dx_uy(1:n1,1:n2),dy_ux(1:n1,1:n2),n1,n2)
    endif
  endif
!
!
!! ------------------------------------------------------------------
!! step iv : find the jacobean in real space 
!! --------------------
!!---------
!$OMP PARALLEL PRIVATE(i1,i2,uu,ww,wXu) SHARED(ukx,uky,omega)
!$OMP DO
  do i2=1,n2
    do i1=1,n1
      uu=0.;ww=0.;wXu=0.
      uu(1)=ukx(i1,i2); uu(2)=uky(i1,i2)
      ww(3) = omega(i1,i2)
      wXu = cross(ww,uu)
! Contribution from immersed boundary and also from the force are added here. 
! At present, we are using the second-order runge-kutta scheme for time stepping
! where the immersed boundary is evolved using an Euler scheme at the end of the
! final step. But, of course, the forces from the boundary points are added in
! both the sub-steps. 
      if (limmersed_boundary) then
!        if (  (i1.eq.138).and.(i2.eq.64)  ) then
!        write(*,*)'DM before',wXu(1:2)
        call add_boundary_force(wXu(1:2),i1,i2)
!        write(*,*)'DM after',wXu(1:2)
!         endif
      endif
      if(lrspace_forcing) call add_rspace_force(wXu(1:2),i1,i2)
      ukx(i1,i2) = wXu(1); uky(i1,i2) = wXu(2)
    enddo
  enddo
!$OMP END DO
!$OMP END PARALLEL
!! ------------Forward Fourier Transform ----------
  call fft_2d_real2complex(ukx)
  call fft_2d_real2complex(uky)
!! -----------------------------------------------
end subroutine eval_jac
